#include <iostream>
#include <vector>
#include "mpreal.h"

using namespace mpfr;

class LinearEquationSolver
{
public:
    LinearEquationSolver(int n)
    {
        this->dims = n;
        A.resize(n*(n+1));

        row_ptrs.resize(dims);
        for (int row = 0; row < dims; row++)
        {
            row_ptrs[row] = &(A[row*(n+1)]);
        }
    }

    void setCoefficient(int row, int col, mpreal value)
    {
        row_ptrs[row][col] = value;
    }

    void setCoefficient(int row, int col, double value)
    {
        row_ptrs[row][col] = (mpreal)(value);
    }

    void set_rhs(int row, double value)
    {
        row_ptrs[row][this->dims] = (mpreal)(value);
    }

    void set_rhs(int row, mpreal value)
    {
        row_ptrs[row][this->dims] = value;
    }

    std::vector<mpreal> solve()
    {
        int N=this->dims;

        std::vector<mpreal> xs;
        xs.resize(N);

        for (int row = 0; row < N; row++)
        {
            int maxRow = row;
            for (int k = row + 1; k < N; k++)
            {
                if (abs(row_ptrs[k][row]) > abs( row_ptrs[maxRow][row]))
                {
                    maxRow = k;
                }
            }

            mpreal* temp  = row_ptrs[row];
            row_ptrs[row] = row_ptrs[maxRow];
            row_ptrs[maxRow]= temp;

            for (int k = row + 1; k < N; k++)
            {
                mpreal factor = row_ptrs[k][row] / row_ptrs[row][row];
                for (int j = row; j <= N; j++)
                {
                    row_ptrs[k][j] -= factor * row_ptrs[row][j];
                }
            }
        }

        for (int i = N - 1; i >= 0; i--)
        {
            xs[i] = row_ptrs[i][N] / row_ptrs[i][i];
            for (int k = i - 1; k >= 0; k--)
            {
                row_ptrs[k][N] -= row_ptrs[k][i] * xs[i];
            }
        }
        return xs;
    }

    int getDims() const
    {
        return dims;
    }

private:
    int dims;
    std::vector<mpreal> A;
    std::vector<mpreal*> row_ptrs;

};

